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Abstract 

In 1918 Brouwer considered stability of a heavy particle in a rotating vessel. This was the first demon- 
stration of a rotating saddle trap which is a mechanical analogue for quadrupole particle traps of Penning 
and Paul. We revisit this pioneering work in order to uncover its intriguing connections with classical rotor 
dynamics and fluid dynamics, stability theory of Hamiltonian and non-conservative systems as well as with 
the modern works on crystal optics and atomic physics. In particular, we find that the boundary of the 
stability domain of the undamped Brouwer's problem possesses the Swallowtail singularity corresponding to 
the quadruple zero eigenvalue. In the presence of dissipative and non-conservative positional forces there is 
a couple of Whitney umbrellas on the boundary of the asymptotic stability domain. The handles of the um- 
brellas form a set where all eigenvalues of the system are pure imaginary despite the presence of dissipative 
and non-conservative positional forces. 

Keywords: Rotating saddle, Quadrupole trap, Crystal optics. Rotor dynamics. Stokes waves, Benjamin- Feir 
instability. Dissipation-induced instabilities. Swallowtail, Whitney umbrella. Exceptional point 



1 Introduction 

In 1918 Luitzen Egbcrtus (Bertus) Jan Brouwer (1881-1966) — a founder of modern topology who established, 
for example, the topological invariance of dimension and the fixpoint theorem — considered a problem on 
stability of a heavy particle in a rotating vessel [T]. The work written in Dutch was almost forgotten until in 
1976 Bottema drew attention to it [2]. However, even those later publications did not manage to popularize 
substantially the Brouwer's work, which is remarkable in many respects. 
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Figure 1: A heavy particle in a rotating vessel [I]. 



In Brouwer's setting, a heavy particle of unit mass is moving on a surface that has a horizontal tangential 
plane at some point O and rotates with a constant angular velocity D, around a vertical axis through O, Fig. [T] 
A rectangular coordinate system {x, y, z) with the origin at O and the [x, z)- and {y, z)-plane of which coincide 
with the principal normal sections of the surface at O, rotates with the surface's angular velocity around the 
z-axis in a counter-clockwise direction, Fig. [1] 

Taking into account that a potential function of the particle is = gz and a kinetic energy is T = ^[(i — 
y^Y + (y + xVtY + z^], Brouwer derives the non-linear equations of motion both in the frictionless case and in 
the case of Coulomb friction between the particle and the surface. The frictionless equations linearized around 
the equilibrium position O, have the form [I] 



where dot indicates time differentiation, fci^2 = 9/fi,2i and ri^2 are the radii of curvature of the intersection 
of the surface with the {x,y)- and (x, z)-plane, respectively, ^ ^2. A merit of Brouwer is that he not only 
allowed the radii of curvature to be non-equal, which was already the state-of-the-art of his time, but he also 
did not impose restrictions on the sign of ki and ^2 and considered both concave (Fig. [Ija)) and saddle-like 
(Fig. mb)) surfaces, ahead of his contemporaries. 

Analyzing the characteristic equation of the system ([T]) 



X - 2Vly + (fci - r2^) x 








0, 



(1) 



A'' + (fci + fc2 + 2n^)\^ + (fci - 02)(fc2 - ^f) ^ 0, 



(2) 
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Figure 2: (a) Contours of the Whitney umbreha and two parabohc cyhnders ^ in the (fci, fc2)-plane for 
17 = 0, 0.3, 0.5, 0.7 and (b) the domain of the marginal stabihty in the same plane for il — 0.7. 

Brouwer writes down the conditions for marginal stability of the particle 

{ki - n'^){k2 - n"^) > 0, 

fci + fca + > 0, 
{ki - k2)^ + 8n\ki + k2) > 0. (3) 

It is instructive to plot the stability domain ([3]) in the (fci, ^2, ri)-space. We restrict our study by the 
non- negative values of fl because the stability domain is symmetric with respect to the plane = 0. 

First we note that the parabolic cylinders ki — ft^ and ^2 — intersect along the parabola that lies in 
the plane fci — k2- This parabola smoothly touches the parabolic cylinder fci + fc2 + 2il^ = as well as the 
Whitney umbrella fci + fc2 = — (fci — k2)'^ / {8fP) at the origin. The Whitney umbrella smoothly touches both 
of the parabolic cylinders fci = Q,^ and fc2 = il^ along spacial curves that project into the lines fc2 = — 3fci and 
fci = — 3fc2 in the (fci, fc2)-plane, see Fig.[2ja). 

In Fig. [2][a) , the contours of the Whitney umbrella and of the two parabolic cylinders are shown for various 
values of il. One sees that the contours of the Whitney umbrella are parabolas that shrink into a fragment of a 
single line fci = fc2 that lies in the third quadrant of the (fci, fc2)-plane. The contours of the parabolic cylinders 
are straight lines parallel to the coordinate axes. The parabolas touch the straight lines, the latter intersect 
each other along the line fci = fc2. Therefore, the stability domains have the form of a curvilinear triangle ABC 
with cuspidal point singularities at the two corners A and B and with the transversal intersection singularity 
at the third one marked as the point C in Fig. ^h) . At the singular point C the triangular stability domain 
{stability for large angular velocities according to Brouwer) is joint with the infinite stability domain given by 
the inequalities fci > Q,^, k2 > 51^ {stability for small angular velocities according to Brouwer). The latter part of 
the stability diagram was found already in 1869 by Rankine [3] whereas the triangle ABC was reconstructed by 
the joint efforts of Foppl (1895) 4 , von Karman (1910) [51, Prandtl (1918) [6], Brouwer (1918) [1,, and Jeffcott 
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Figure 3: The Swallowtail sing ularity (0*, TT) at the origin on the stability boundary of the Brouwer's problem 
on a heavy particle in a rotating vessel illustrating peculiarities of gyroscopic stabilization. 

(1919) 7 . Note that similar planar stability domains appear in a recent study of levitation of a rotating body 
carrying a point electrical charge in the field of a fixed point charge of the same sign i8J , see also |9] . 

A remarkable result of Brouwer is that for kik2 < 0, i.e. when the surface is a saddle, the particle can be 
stabilized by rotation of the surface. Brouwer concludes that 'as long as in O the principal curvature, concave 
in an upward direction, is less than three times the one concave in a downward direction, there are rotation 
velocities for which the motion of the particle on the rotating saddle yields formal stability' T . Cuspidal points 
A and B on the stability boundary clearly seen in Fig. EJb) evidence that such a stabilization occurs for a rather 
non-trivial choice of the main curvatures of the rotating saddle. 

We notice however, that while Brouwer [1^ neither plotted the stability domain nor discussed its geometry, 
Bottema 2^ improved Brouwer's result and even produced a planar stability diagram which however does not 
detect the singularities because of the unsatisfactory choice of parameters. Therefore, we substantially refine 
the analysis performed by Brouwer and Bottema by finding the singularities of the planar stability boundary for 
the system ([1]). Moreover, we discover that in the {ki, fc2, ri)-space the full stability domain given by inequalities 
^ is bounded by a surface that has the Swallowtail singularity at the origin, see Fig. [3l Bold lines in Fig. [3] 
indicate two cuspidal edges (A, B) and one regular intersection (C) on the stability boundary. Thus, the part of 
the stability domain that corresponds to stability for large angular velocities is inside the spike of the singular 
surface whereas the part of the stability domain that corresponds to stability for small angular velocities occupies 
the remaining bigger portion of it that has only a regular edge singularity. The tangent cone to the stability 
domain at the tip of the spike consists of the semi-axis 51 > 0, Fig. [3] Similar structure of the stability domain 
was found in the model of a two- link rotating shaft with four degrees of freedom [Qj [10] . 

At the origin, the conditions fci = fc2 = and f2 = yield quadruple zero eigenvalue A = with the Jordan 
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block. Galin [TT had proven that this spectral degeneracy (0^) corresponds to the Swallowtail singularity of 
the stability boundary of a Hamiltonian system and has codimension 3. The Swallowtail singularity generically 
appears in the Brouwer's problem that is described by the three parameter Hamiltonian system ([1}. We note 
that in the works [9l I10| the Swallowtail singularity had been found on the stability boundary of a two-link 
rotating shaft that has twice as many degrees of freedom. The Brouwer's problem on a heavy particle in a 
rotating vessel appears to be the simplest physically meaningful system that possesses this singularity and as 
we will show further in the text, this fundamental paradigmatic model has deep connections to many important 
effects of classical and modern physics. 

2 Manifestations in physics 

The singular surface bounding the stability domain of the Brouwer's heavy particle in a rotating vessel perfectly 
symbolizes the variety of applications to which it is connected. Almost literally at every corner of this surface 
lives a physical phenomenon. 

2.1 Helical quadrupole magnetic focussing systems 

Already in 1914, in his correspondence to Georg Hamel, Brouwer discussed the possibility of experimental 
verification of stabilization of a heavy ball on a rotating saddle [1] . Nowadays such mechanical demonstrations 
are available as teaching laboratory experiments [12 . However, the prospects for the first non-trivial physical 
application of this effect had arisen as early as 1936 when Penning proposed using the quadrupole electric field 
to confine the charged particles [131 [H] . 

In the RF-electric-quadrupole Paul trap invented in 1953 [15l [16], a saddle-shaped field created to trap 
a charged ion is not rotating about the ion in the center. The Paul-trap potential can only 'flap' the field 
up and down, which yields two decoupled Mathieu equations describing the motions of a single ion in the 
trap. Nevertheless, comparing the rotating saddle trap and the Paul trap, Shapiro |171 118] and Thompson 
et al. [12] demonstrated that the former mimics most of the characteristics of the Paul trap such as stability 
and instability regions, micromotion and secular oscillation frequency 12^. Despite the Brouwer's problem is 
not only a mechanical analogy to the Paul trap, the equations of the two models are substantially different. 
Brouwer's system is autonomous, which greatly simplifies its solution. An example of a natural appearance 
of the Brouwer's equations comes from accelerator physics, where they originate in a theory of focussing by a 
helical magnetic quadrupole lens [19]. 

Indeed, according to Pearce jl9j the linearized equations of motion in the laboratory (x, y, z)-frame of a 
particle of momentum p and charge e in the helical quadrupole magnetic field that rotates with distance along 
the z-axis, completing one rotation in an axial distance A, in the assumption that x,y X, are 

X — a(x cos z -\- y sin z) , 

y = a(a;sinz — j/cosz), (4) 
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where dot stands for the derivative d/dz and 



Yk^ , 2 eG 



and c is the speed of light. The downstream distance z is measured in units of A/(47r). At z = 0, the components 
of the magnetic field in the laboratory frame are = Gy, By — Gx, i?^ = as for the conventional quadrupole. 
By coordinate transformation to a frame {X, Y, z), which rotates with the mechanical rotation of the helix: 

z . z 
X — xcos — hysin— , 
2 ^ 2' 

z z 

Y = — a; sin — + y cos — , (6) 
we arrive at the autonomous system of equations that is equivalent to ^ 

v + Jv+(K+(J/2)2)v = 0, (7) 

where 

We note that (l7|) is a particular case of the Brouwer's rotating vessel equations ([T]) with fi = 1/2 and k2 = 
— fci — a that corresponds to a rotating saddle. 

2.2 Propagation of polarized light within the cholesteric liquid crystal 

Following Marathay [20] we consider the state of polarization of a plane wave as it propagates along the helical 
axis within the nonmagnetic cholesteric liquid crystal. Assume that z-axis of a right-handed laboratory frame 
is parallel to the helical axis. The molecular planes of the structure are parallel to the {x, y)-plane of the 
laboratory frame. The principal axes of each molecular plane are gradually rotated as one proceeds along z- 
direction. Define the pitch p as the minimum distance between the two planes whose principal axes are parallel 
and introduce the parameter a — 2tt/p. Restricting ourselves by the plane waves that propagate normally to 
the crystal planes, i.e. along the helical axis, we find from the Maxwell's equations that 

where E ~ E(z,t) is the Jones vector [20]. The 2x2 matrix e{z) describes the dielectric properties of the layer 
at z in the laboratory frame and for a right-handed structure in terms of the principal dielectric constants eli 
and it can be expressed as 

e(z) = R(az)KR-i(az), (10) 
where K = diag (en, £§2) and R is the rotation matrix 

. cosfaz) — sm(az) . 

R= I . I • (11) 

sin(az) cos(az) 
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Looking for a time-harmonic solution of Eq. ^ of the form E(z, t) — F {z)exp{—iujt) and transforming to a 
space-rotating coordinate system by 

f(z) = R(-az)F(z), (12) 

we find that in the space rotating coordinate system the dynamics of the state of polarization f is governed by 
the autonomous equation 

f + 2aJf +(fc2K + (aJ)2)f = 0, (13) 

where fc^ ~ uj^ jc? and the dot stands for the derivative djdz. We again obtained the Brouwer's equations ([T]), 
with f7 = a, fci = fc^e"]^ and k2 = ^^£22- 

2.3 Benjamin-Feir instability of traveling waves 

The modulations of Stokes waves in deep water are described by the Non-Linear Schrodinger equation (NLS) 

lAt-l-aA,, -f 7|^|2^ = 0, (14) 

where i — \/— T, A is the envelope of the wave carrier, a and 7 are positive real numbers and the modulations 
are restricted to one space dimension a;, for details see recent work by Bridges and Dias [21]. 

Assuming A = u\ -I- 11/2, linearizing the non-linear partial differential equation (|14|) about the basic traveling 
wave solution with the amplitude Uo ~ (m", ^2)1 spacial wavenumber k and the frequency uj and then expanding 
the periodic in x solution with the period cr of the linearized problem into Fourier series, we find that the 
(T-dependent modes decouple into four-dimensional subspaces for each harmonic with the number n, so that for 
n = 1 we get 

Jv + 2Q!fc(TJw — acr^v + 27UoU(^v = 0, 

Jw — 2akaJ'v — aa^w + 27UoU^w = 0, (15) 

where dot indicates time differentiation, the matrix J is defined in (jS)), the dyad uqUj^ is a 2 x 2 symmetric 
matrix and — ak^ — ||uo|p. The NLS model ([T4l) has a Benjamin-Feir instability: considering all other 
parameters fixed, with the increase of the wave amplitude Uq at some threshold value two eigenvalues of the 
linear stability problem (jlSp experience the Krein collision [22^ , and these two modes have negative and positive 
energy, respectively [21]. 

Differentiating the first of the equations ()15|) once, substituting into the result the expression for Jw that 
follows from the second of the equations (jlSp. then expressing w from the first of the equations ()15|) and finally 
multiplying the result by J from the left we obtain the equation that describes the time evolution of v 

v-2JBv + (4a^fcV^ + (JB)2)v = 0, (16) 

where B — 27UoU^ — aa^l, and I is the unit matrix. The vector w satisfies the same equation. The matrix JB 
is non-symmetric and thus can be decomposed into a symmetric and skew-symmetric parts 

JB = (zJ + 7D, (17) 
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where q — a^a — 7||uo|p and D — UquJ^J — JuqU;^. 

With this notation we reduce the equation ([TO]) to the standard for the gyroscopic systems form 



V + 2gJv + 27DV + (P + (gJ)2)v = 0, 



(18) 



where P = {'ia^ k'^ a'^ + 72||uo||'')I. 

Equation (fT8|) is the damped Brouwer's equation ([T]) with VL = q and fci = fc2 = 4Q;^A:^cr^ + 7^||uo||^. 
The 'damping' matrix D is symmetric and traceless, i.e. it is indefinite with the eigenvalues /ii^2 = i||uo|p. 
Therefore we can give a new interpretation to the Benjamin- Feir instabihty and treat it as the destabilization 
of a gyroscopic system by the indefinite damping [23 . Remarkably, this destabihzation mechanism is known in 
rotor dynamics, where the indefinite damping matrix results from the falling friction characteristic that typically 
occurs in the problems of acoustics of friction such as the squealing brake or the singing wine glass [531 IMl US] . 
In [25j destabilization of the Brouwer-type equations by the indefinite damping was considered in detail. 

2.4 Stability of deformable rotors 

With the equal positive coefficients fci = fc2 > 0, the Brouwer's equations ([T]) coincide with that of the idealized 
model of the classical rotor dynamics problem of shaft whirl by Foppl (1895) |4j, von Karman (1910) and 
Jeffcott (1919) \Jj, written in the rotating (x, ?;)-frame. 

Indeed, a deformable shaft carrying, e.g., a turbine wheel, and rotating with the angular speed about its 
axis of symmetry can be modeled as a planar oscillator on a rotating plate [26j . i.e. as a unit mass point that is 
suspended symmetrically by massless springs with the effective stiffness coefficients ki — k2 — k from the frame 
that rotates with the angular velocity [271 [28], see Fig.[4l^a). 

For the symmetrical Foppl- von Karman- Jeffcott rotor the transverse bending modes occur in pairs with equal 
natural frequencies a;„ — \fk for vibrations in orthogonal diametral planes. When the vibrations of such a pair 
are combined with equal amplitudes and a quarter period phase difference, the result is a circularly polarized 
vibration — a clockwise- or a counter-clockwise circular whirling motion with the whirl rate |27[ I28| . This 
phenomenon is related to the Doppler splitting of the doublet modes into the forward and backward traveling 
waves, which propagate along the circumferential direction of a rotating elastic solid of revolution discovered 
by Bryan in 1890 [29] . 




According to the Brouwer's stability conditions ([3|) the symmetrical Foppl-von Karman-Jeffcott rotor is 
marginally stable at any speed 51. Indeed, the part of the plane k\ — k^ corresponding to k\^2 > belongs to 
the domain of marginal stability, as is seen in the stability diagrams of Fig. ^h) and Fig. [3| 

The symmetry of the Foppl-von Karman-Jeffcott rotor is, however, a latent source of its instability. Destabi- 
lization can be caused already by constraining the mass point to vibrate along a rotating diameter. For example, 
if a guide rail is installed along the a;-axis so that y is constrained to vanish identically, equation ([1]) reduces to 
the 1869 model by Rankine [3] 




x+iki~ n^)x = 0, 



(19) 
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Figure 4: (a) A Prandtl-Smith model of whirling shaft [371 US] : A point mass m = 1 suspended by the springs 
from the frame that rotates with the angular velocity il. (b) Whirling frequencies of the Prandtl-Smith rotor 
with fci 1 as a function of fl and e form a singular surface with four conical points in the plane e = 0. (c) 
Growth rates of the Prandtl-Smith rotor showing instability due to stiffness detuning near the critical speeds 

n = ±1. 

which predicts unbounded growth for x as soon as > ki, i.e. as soon as the centrifugal field overpowers 
the elastic field [13 [2B]. In the diagram of Fig. [IJb) Rankine's instability threshold (critical rotating speed) 
bounds the infinite stability domain with the corner at the point C showing that the constraint could also be 
achieved by tending one of the stiffness coefficients to infinity. Note that recently Baillieul and Levi motivated by 
stabilization of satellites with flexible parts, studied the dynamical effects of imposing constraints on the relative 
motions of component parts in a rotating mechanical system or structure and in particular in the Brouwer's 
equations i30l . 

Although after successful demonstration by De Laval of a well-balanced gas turbine running stably at 
supercritical, i.e. unstable by Rankine, speeds in 1889, the model of Rankine [3] has been recognized as 
inadequate |271 128] . we see that its instability threshold is just a part of the singular stability boundary of the 
Brouwer's equations shown in Fig. [H 

Soon after the success of the symmetric Foppl-von Karman-Jeffcott model of De Laval's rotor, Prandtl 
(1918) [6 noted that if the elastic or inertia properties of a rotor are not symmetric about the axis of rotation 
there may be speed ranges where the rotor whirl is unstable. He pointed out the analogy between a rotor 
whose shaft had unequal principal stiffness coefficients (fci ^ and a pendulum mounted asymmetrically on a 
turntable. The pendulum was unstable when the turntable rotation rate 57 was between the pendulum natural 
frequencies in the soft and stiff directions. Similar instability due to inertia asymmetry was discussed by Smith 
in 1933 [31]. For positive coefficients fci 7^ fc2, the Brouwer's equations ([T]) is exactly the undamped model by 
Prandtl and Smith. 

Let us assume in the equation ^ that fci = 1 + e and fc2 = 1 — £ [271 [53] . Calculating the imaginary parts 
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of the roots of the characteristic equation ^ we find the imphcit equation for the whirhng frequencies in the 
Prandtl-Smith model 

((ImA)2 - 1 - 172)2 _ 4f^2 ^ ^20) 

In the e, ImA)-space the whirhng frequencies he on a singular surface that has four conical points, Fig. |4ljb). 
The slice of this surface by the plane e — shows four straight lines that intersect each other at the conical 
points at = and at the critical speeds = ±1. These eigencurves constitute the Campbell diagram [32] 
of the ideal Foppl-von Karman-Jeffcott rotor. Stiffness modification with the variation of s corresponds to the 
eigencurves in a slice of the eigenvalue surface by the plane that departs of the conical singularities by a distance 
£ 7^ 0. In the vicinity of = the eigencurves demonstrate avoided crossings while near the critical speeds 
fl — ±1 the eigenvalues are real with zero frequency. As is seen in Fig. SJc), the real eigenvalues (growth rates) 
form a singular eigenvalue surface with two conical points at the critical speeds 

((ReA)2 + 1 + VL^f - 4f|2 ^ ^21) 

The slices of the surface ([21]) by the planes e 7^ in the vicinity of the apexes of the cones are closed loops 
known as the 'bubbles of instability' |33j. Hence, for small £, the Prandtl-Smith rotor is unstable at the speeds 
in the interval 

< \^\ < VTTs. (22) 

The amount of kinetic energy available in the rotor is usually orders of magnitude greater than the de- 
formational energy which any internal mode can absorb. Therefore, already small deviations from the ideal 
conditions yield coupling between modes that by transferring even a tiny fraction of the rotational energy can 
initiate failure in the vibratory mode [371 [2E] ■ 

The conical singularities on the surfaces of frequencies and growth rates arise as unfoldings of double semi- 
simple eigenvalues of the Foppl-von Karman-Jeffcott rotor that exist at £7 = and at the critical speeds. In 
1986, MacKay [33^ studied unfoldings of the semi-simple double pure imaginary and zero eigenvalues in general 
Hamiltonian systems and demonstrated that the orientation of the cones depends on the Krein signature |22j 
of the doublets. At f2 = the doublets have definite Krein signature and according to MacKay the frequency 
cones are oriented vertically in the (il, e, ImA)-space that yields avoided crossings in their slices and thus stability 
[231 133] ■ If the doublet has a mixed Krein signature then there appears a differently oriented eigenvalue cone for 
the real parts of the perturbed eigenvalues that possesses bubbles of instability in its cross-sections [24l [33l [34] . 

According to ()22|) the whirling instability onset is no less than the critical speed for the mode in question. 
This yields a traditional for rotor dynamics suggestion that 'although Rankine was wrong in thinking that 
supercritical operation was always unstable, he was partially correct in concluding that dynamic instability 
would not occur at subcritical speeds' [28] . 

Indeed, for perturbations that preserve the Hamiltonian structure of the equations of motion of the gy- 
roscopic system ([T]), the instabilities in the subcritical speed range are prohibited by the definiteness of the 
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Krein signature of the doublets. However, the non-Hamiltonian perturbations caused by dissipative and non- 
conservative positional forces may yield dynamic instabilities at the small rotational speeds in the subcritical 
speed range as it happens, for example, in the problems of acoustics of friction |24) . Subcritical dynamic insta- 
bilities have many unusual properties partially uncovered in recent works [24i I25i 134) . In the next section, we 
report new effects in the subcritical speed range on the example of the near-Hamiltonian Brouwer's problem. 

3 Brouwer's system in the presence of dissipative and non- conservative 
positional forces 

In 1976 Bottema [5] extended Brouwer's original setting by taking into account the internal and external damping 
with the coefficients i5 and i/, respectively 



where z — {x, y)'^ , D — diag {S + 1^,6 + i'), G = J, K = diag (fci, fc2), and N = QJ. 

Note, however, that the equation ((23)) appears first in a 1933 work by Smith [31] as a model of a rotor 
carried by a flexible shaft in flexible bearings, where stationary (in the laboratory frame) damping coefficient v 
represents the effect of damping in bearing supports while rotating damping coefficient 6 represents the effect 
of damping in the shaft. These two types of damping were introduced in 1924 by Kimball |35| in order to 
explain destabilization due to damping within a rotor, studied also by Kapitsa 36^. In a more general model 
of a rotating shaft by Shieh and Masur [37] the damping matrix is allowed to be D = diag (i5i, (^2) while the 
matrix of the non-conservative positional forces is simply N = J. In recent works the matrices D 

and K were assumed to be arbitrary symmetric while N = J with v — ^'(f2) as an arbitrary smooth function of 
the rotation speed. 

The system is a particular case of a general non-conservative system with two degrees of freedom 



where A and B are real non-symmetric matrices. Its characteristic equation is given by the Leverrier-Barnett 
algorithm |38j 



Z + BZ + 2r2Gz + (K + (f7G)2)z + j/Nz = 0, 



(23) 



z -I- Bz -f Az = 



(24) 



A'* + trBA^ + (trA + dct B)A^ + (tr AtrB - trAB)A + det A = 0. 



(25) 



The equation (l25t becomes biquadratic when the following two conditions are fulfilled 



trB = 0, trAB = 0. 



(26) 



Since trAB = trBA, these conditions are reduced to 



tiBs = 0, trA^Bs + trAaBa = 



(27) 
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Figure 5: Imaginary and real parts of the roots of the characteristic equation (|28p as a function of the damping 
coefficient under the constraints ((29|) for ki — 1, ft — 0.03 and v — 0.03. 

where subscripts s and a denote symmetric and skew-symmetric components of the matrices A and B. The 
conditions (|27|) are satisfied in particular in the case of the gyroscopic system when = and Aa = 0, and in 
the case of the circulatory system when B^ = and Ba = 0. 

The systems, different from circulatory and gyroscopic that nevertheless satisfy conditions ((27)) should have 
some unusual properties such as an indefinite damping matrix with zero trace. Although at a first glance such a 
physical system looks exotic, the non-trivial example is given already by the classical Benjamin-Feir instability 
of the Stokes waves discussed in the previous section. Indeed, in the equation psp the symmetric matrix D is 
traceless and the matrix of positional forces P + (qJ)^ is proportional to the unit matrix which yields fulfillment 
of the conditions (P7)) . As a consequence, the eigenvalues of the system (IT51) are symmetric with respect to 
the imaginary axis in the complex plane which implies marginal stability if they are pure imaginary. The loss 
of stability of the system (jlSp with the indefinite damping happens through the collision of two imaginary 
eigenvalues into a double eigenvalue with the Jordan block and subsequent splitting of it into two complex ones 
with positive and negative real parts (Krein collision). This is not surprising in light of that the system (IT51) 
was shown to be Hamiltonian in |21) . 

The model of Shieh and Masur [37] provides a more non-trivial example where the constraints ([27]) play a 
role. Indeed, its characteristic equation has the form 

X^ + {6i+62)X^ + {6i62 + ki+k2+2n'^) + (fci^a + ^1^2 + 40i/ - (<5i -t- Jz)!^^) X + i^l"^ -ki){n^ - k2) + iy^ ^ 0. 

(28) 

The equation is biquadratic in case when 

5i+S2=0, ^ = -^, (29) 

Ol 
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with K ^ k2 — ki. If fci > r2^ and v > then all the roots of the equation (j28p are pure imaginary when 

In Fig. [5] the pure imaginary eigenvalues are shown by black lines as functions of the damping parameter Si. At 

^' = ^' .- + {k,-n^r " = + " "fcT^ ^^^^ 

there exist a double zero eigenvalue with the Jordan block (0^), see Fig. [5] and Fig. EJa). In the interval 
Q < 5i < 5d there exist one positive and one negative real eigenvalue. In Fig. [5] the eigenvalues with non-zero 
real parts are shown in red. In the ((5i,(52, K)-space the exceptional points (EPs) 

{-2Vl,2Q.,2v), (217,-217,-21/) (32) 

correspond to the double pure imaginary eigenvalues with the Jordan block 



A_20 = ±^v/fcl - + jy, = ±i^/ki - - (33) 

for 5i = —217 and 5i ~ 217, respectively. 

We see in Fig. [5] that changing the damping parameter 5i we migrate from the marginal stability domain 
to that of flutter instability by means of the Krein collision of the two simple pure imaginary eigenvalues as 
it happens in gyroscopic or circulatory systems without dissipation. It is remarkable that such a behavior of 
eigenvalues is observed in the gyroscopic system in the presence of dissipative and non-conservative positional 
forces. The existence of pure imaginary spectrum in the non-conservative dissipative systems satisfying the 
constraints (|27l) reminds the similar phenomenon in PT-symmetric non-Hermitian systems that can possess 
pure real spectrum [39]. 

Let us now establish how in the ((5i, (52, '^)-space the domain of marginal stability given by the expressions 
(Uni) and ([50]) is connected to the domain of asymptotic stability of the non-constrained equation (pS)) . Writing 
the Lienard and Chipart [40] conditions for asymptotic stability of the polynomial (|28p we find 



pi := 81+82 > 0, 
P2 := 8182 + ki + k2 + 2172 > 0, 

P4 -.^ {n^ - ki)(n^ - k2) + > 0, 

H3 := (<5i + S2)iSi82 + ki+ k2 + 2172)(A:i(52 + 5ik2 + 417i/ - {81 + ^2)^^') 
-{81+82? {{n^ -ki){n^ ~k2) + v'') ~{ki82 + 8ik2+ mv- {81+82)^1^ f > 0. (34) 

The surfaces pi — Q and Hj, = are plotted in Fig. \6la) . The former is simply a horizontal plane that passes 
through the point of the double zero eigenvalue (0^) with the coordinates {8d, —8d, k^.) and thus bounds the 
stability domain from below. The surface Ht, = is singular because it has self-intersections along the portions 
of the hyperbolic curves ([29]) selected by the inequalities (l30l) . The curve of self-intersection that corresponds to 
K > ends up at the EP with the double pure imaginary eigenvalue A_20- Another curve of self- intersection has 
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Figure 6: Stability domain of the rotating shaft by Shieh and Masur for fci = 1, O = 0.03, v = 0.03. (a) The 
'Viaduct' in the {6i, 62, K)-space and its slices in the {61, (52)-plane with (b) k ~ —0.06 and (c) k — —0.03. 

at its ends the EP with the double pure imaginary eigenvalue X2n and the point of the double zero eigenvalue 
(0^). In Fig. [6l^a) the curves of self-intersection are shown in red and the EP and 0^ are marked by the black 
and white circles, respectively. At the point 0^ the surfaces p4 = and = intersect each other forming 
a trihedral angle singularity of the stability boundary with its edges depicted by red lines in Fig. [Hfa). The 
surface H3 = is symmetric with respect to the plane pi — 0. Thus, a part of it that belongs to the subspace 
Pi > bounds the domain of asymptotic stability. 

At the EPs, the boundary of the asymptotic stability domain has singular points that are locally equivalent 
to the Whitney umbrella singularity [41] |42j |43l |44] |45]. Between the two EPs the surface H3 = has an 
opening around the origin that separates its two sheets. This window allows the flutter instability to exist in 
the vicinity of the origin for small damping coefficients and small separation of the stiffness coefficients k. In 
[34j an eigenvalue surface of the same shape was named the Viaduct. 

In Fig. IHb) a cross-section of the surface — by the horizontal plane that passes through the lower 
exceptional point is shown. The domain in grey indicates the area of asymptotic stability. Its boundary has a 
cuspidal point singularity at the EP. Although the very singular shape of the planar stability domain is typical 
in the vicinity of the EP with the pure imaginary double eigenvalue with the Jordan block 45J, the unusual 
feature is the location of the EP that corresponds to non- vanishing damping coefficients, Fig. Mjo). 

Indeed, in the known studies, the system with vanishing dissipation was either gyroscopic (Hamiltonian) 
as in the works [40l |45l |46l |47l SH |49l [50] or it was circulatory (reversible) as in the works [HI [45l [46l [51] . 
The undamped system in such problems has a domain of marginal stability bounded by the exceptional points 
that correspond to double pure imaginary eigenvalues. This is not the case in the Brouwer's problem with the 
dissipative and non-conservative forces by Shieh and Masur [37]. According to the theorems of Bottema [52] 
and Lakhadanov [53] the undamped gyroscopic system with non-conservative positional forces is generically 
unstable. By examining the slices of the surface H3 = at various values of k one can see that the origin is 
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indeed always unstable, Fig. infb,c). At k = the origin is unstable in the presence of the non-conservative 
positional forces even when the rotation is absent (ft — 0) according to the theorem of Merkin 49 . Contrary to 
the situation known as the Ziegler's destabilization paradox [51 , in the Brouwer-Shieh-Masur model the tending 
of the damping coefhcients to zero along a path in the (52)-pl£ine cannot lead to the set of pure imaginary 
spectrum of the undamped system because in this model such a set corresponds to the non-vanishing damping 
coefficients. 

On the other hand, the lower bound on the difference of the coefficients fci_2 given by the equation (|3ip 
depends on the magnitude of the non-conservative positional force v. li v = 0, the heavy particle on a rotating 
saddle is destabilized by the arbitrary small internal damping, because in this case Kd = —ki+Q? and the 
stability condition k > Kd reduces to K2 > il^, in full agreement with [2]. Nevertheless, with the increase of v 
the critical difference of ki and k2 takes bigger negative values so that for ly^ > il'^(ki — il^) asymptotic stability 
is possible also for the coefficient k2 < while fci was assumed to be positive. 

Therefore, the Brouwer-Shieh-Masur model provides a non-trivial example of a gyroscopic system that can 
have all its eigenvalues pure imaginary in the presence of dissipative and non-conservative forces. The corre- 
sponding set in the three-dimensional space of damping coefficients and the difference between the curvatures 
(stiffness coefficients) is a space curve that ends up at the exceptional points of double pure imaginary eigenval- 
ues with the Jordan block. The set forms singularities of the boundary to the domain of asymptotic stability 
of the Brouwer-Shieh-Masur model such as Whitney umbrellas and dihedral and trihedral angles. The highly 
non-trivial shape of the discovered stability boundary illustrates the peculiarities of destabilization of a rotor in 
the subcritical speed range as well as the stabilization of a heavy particle on a rotating saddle for high rotating 
velocities by non-conservative positional forces in their interplay with the dissipative, gyroscopic and potential 
ones. 

4 Conclusion 

The Brouwer's problem on a heavy particle in a rotating vessel is a paradigmatic model for studying instabilities 
in diverse areas of physics. Brouwer-like systems of equations naturally arise for instance in accelerator physics, 
crystal optics, theory of Stokes waves in deep water, as well as in classical rotor dynamics. Despite a long 
history, the model of Brouwer was not completely understood until now. We established new results on stability 
in the Brouwer's problem that include discovery of the Swallowtail singularity on the stability boundary of 
the undamped Brouwer's equations. This singularity was overlooked both by Brouwer and by Bottema who 
re-visited Brouwer's work in 1976. Moreover, we studied Brouwer's equations in the presence of dissipative and 
non-conservative forces and found a new type of the asymptotic stability domain in the three-dimensional space 
of parameters. It turns out that the stability boundary is a surface with two Whitney umbrella singularities, 
the 'handles' of which correspond to systems with the pure imaginary spectrum despite the presence of damping 
and non-conservative positional forces. We established the necessary conditions for a non-conservative system to 
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have a pure imaginary spectrum and found that the hnearized equations describing the onset of the Benjamin- 
Feir instabiUty satisfy them. As a consequence, the Benjamin-Feir instabihty was interpreted as destabiUzation 
of a gyroscopic system by indefinite damping. 

To summarize, the Brouwer's model describes two common situations: a heavy particle in the narrow concave 
rotating vessel is difiicult to destabilize; that one in the wide rotating dish or on a rotating saddle is difficult to 
stabilize. The former is common in the modern problems of low-speed rotor dynamics related to squealing brakes 
or singing glasses of the glass harmonica. The latter is typical for the problems of particle confinement in atomic 
and accelerator physics as well as for high-speed rotor dynamics. In both situations the destabilization and 
stabilization mechanisms exist that involve dissipative and other non-conservative forces. The very possibility 
to create (in)stability in such circumstances where it conventionally is not expected makes such problems a 
reach source for new development and insight in stability theory. The Brouwer's problem on a heavy particle in 
a rotating vessel sharply highlights some crucial eff'ects that arise in the field of dissipation-induced instabilities. 
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